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A  three-dimensional  parabolic  equation  (3DPE)  that  handles  wide  angles  in  the  vertical,  narrow 
angles  in  the  azimuth,  and  rough  ocean  surfaces  and  bottoms  is  derived.  The  3DPE  is  solved 
numerically  using  the  method  of  alternating  directions.  Surface  roughness  is  accounted  for  by  a 
reflection  coefficient  that  depends  on  grazing  angle.  Calculations  are  presented  to  illustrate  the 
rough  surface  model  and  to  demonstrate  that  azimuthal  diffraction  can  be  important  in  shallow 
water.  The  ability  of  the  3DPE  to  accurately  handle  azimuthal  diffraction  is  demonstrated  with  a 
benchmark  calculation.  Algorithms  for  improving  the  efficiency  of  3DPE  models  are  discussed.  J. 


1.  Introduction 
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The  parabolic  equation  [l]  (PE)  method  is  a  very  useful  model  for  underwater  propagation 
calculations.  For  most  applications  of  the  PE,  the  domain  is  assumed  to  be  cylindrically  symmetric 
with  smooth  boundaries,  and  the  two-dimensional  PE  (2DPE)  is  applied.  Including  the  effects 
of  three-dimensional  variations  and  rough  surfaces  in  the  model  would  allow  one  to  handle  more 
realistic  problems.  The  three-dimensional  PE  (2-5]  (3DPE)  has  been  derived  for  problems  involving 
variations  in  both  range  and  azimuth.  For  some  problems,  accurate  results  can  be  obtained  for 
problems  involving  variations  in  azimuth  by  using  the  2DPE  in  each  direction  of  interest  [3].  For 
such  problems  the  azimuthal  diffraction  term  in  the  3DPE  is  apparently  less  important  than  the 
azimuthal  refraction  term.  A  3DPE  that  handles  wide-angle  propagation  in  both  depth  and  azimuth 
has  been  derived  [4],  However,  existing  numerical  solutions  of  this  3DPE  are  inefficient  due  to 
coupling  between  the  depth  and  azimuth  term.  A  3DPE  that  handles  wide  angles  in  depth  and 
narrow  angles  in  azimuth  can  be  solved  very  efficiently  with  the  method  of  alternating  directions 
[5],  Propagation  in  the  ocean  can  be  significantly  affected  by  rough  boundaries.  [6]  An  efficient 
approach  for  modeling  rough  ocean  surfaces  and  bottoms  is  to  assume  that  the  net  result  of  the 
roughness  is  an  energy  loss  that  depends  on  grazing  angle  [7,8], 


In  this  paper,  a  3DPE  is  derived  that  handles  wide  angles  in  depth  and  narrow  angles  in  azimuth 
and  includes  rough  boundaries.  Rough  ocean  surfaces  and  bottoms  are  handled  approximately  by 
assuming  a  reflection  coefficient  at  the  ocean  surface  that  depends  on  grazing  angle  and  transforming 
this  into  a  homogeneous  boundary  condition  that  is  easily  incorporated  into  the  numerical  solution 
i  of  the  3PDE.  The  reflection  coefficient  is  allowed  to  vary  in  both  range  and  azimuth.  Calculations 

are  presented  to  illustrate  the  rough  surface  PE  and  to  demonstrate  that  the  2DPE  can  not  handle 
some  shallow-water  propagation  problems.  Efficient  algorithms  for  the  numerical  solution  of  the 
3DPE  are  discussed  and  illustrated  with  examples.  The  ability  of  the  3DPE  to  accurately  handle 
azimuthal  diffraction  is  demonstrated  with  a  benchmark  calculation. 


STATEMEftt  A 

Approved  lor  pottle  lohaeef 
Dtofrt  button  UsJbnMod 


07 


099 


248 


*  A 

M.D.  Collins  and  SA.  Chin-Bing 


2.  The  three-dimensional  parabolic  equation 


We  work  in  cylindrical  coordinates  with  2  being  the  depth  below  the  ocean  surface,  9  being 
the  azimuth  angle,  and  r  being  the  horizontal  distance  (range)  from  a  time-harmonic  point  source 
of  circular  frequency  u>.  For  now,  we  assume  that  the  complex  wavenumber  K  =  k  +  itj/3  |  k  |  and 
the  density  p  depend  only  on  2,  where  fc  =  u/c,  /?  is  the  attenuation  in  decibels  per  wavelength,  c 
is  the  sound  speed,  and  rj  =  (40irlog10e)_1.  We  define  the  reference  sound  speed  Cq  and  reference 
wavenumber  ko  -  w/c0.  Cylindrical  spreading  is  handled  by  removing  the  factor  r~i  from  the 
complex  pressure  P. 

For  kr  >>  1,P  is  assumed  to  satisfy  the  farfield  equation 


_  1  dpap  a^p  j_£p  2 

dz 2  ~pdzdz  +  ar3  +  r*dO*  + 


(2.1) 


Since  r  3  may  be  assumed  to  commute  with  d/dr  in  the  farfield,  we  obtain  the  following  factor¬ 
ization  for  the  outgoing  solution 
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We  define  the  operators 


x  =  fc0-3  (a-3  -  *03  +  ^ 


1  dp  d  \ 

pdzdz  J 


(2.2) 


(2.3) 


Y  = 


1  a 3 

k* r3  as* 


(2.4) 


and  assume  that  |  Y P  |<<|  X P  |.  For  the  square  root  in  Eq.  (2.2),  we  use  the  approximation 


_ _  o  y  1 

VT+X+T  =  1  +  —  +  -Y  +  0(X\XY) 

and  remove  the  plane  wave  factor  exp(tfc„r)  from  P  to  obtain  the  3DPE 
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Figure  1:  The  modal  radiation  patterns  for  Pa  at  r  =  300m  (dashed  curves)  and  for  P  (solid  curves) 
at  (a)  r  =  150m  and  (b)  r  =  300m. 

Equation  (2.6),  which  reduces  to  the  wide-angle  2DPE  [9]  for  two-dimensional  problems,  is  a  valid 
leading-order  solution  to  problems  in  which  K  and  p  depend  on  r  and  6  as  a  perturbation.  An 
analogous  approach  based  on  Eq.  (2.5)  has  been  used  to  derive  a  wide-angle  time-domain  PE  for 
shallow  water.  [10] 

We  solve  the  3DPE  numerically  with  the  method  of  alternating  directions,  [11]  which  requires 
numerical  methods  for  each  of  the  following 


ep  _  **■  (gi  -  +  a?  -  Ilf £}) 
*  =  “S  + 

dP  _  i  d'P 
dr  2 k^r2  dB2 


(2.7) 

(2.8) 


Equation  (2.7)  is  the  wide-angle  2DPE  and  can  be  solved  with  the  numerical  methods  described  in 
Ref.  10.  Equation  (2.8)  can  be  solved  with  centered  differences  in  6  with  Crank-Nicolson  integration 
in  r.  The  matrices  involved  have  three  diagonals  and  entries  in  the  upper  right  and  lower  left  corners 
for  the  continuity  condition 


P  |p=0  =  P  |0=2>r  ■ 


(2.9) 


□ 

□ 


To  show  that  the  3DPE  accurately  handles  azimuthal  diffraction  (i.e.  the  effect  of  the  0-term 
in  the  3DPE),  we  consider  a  dipole  source  with  a  horizontal  polar  axis.  In  an  ocean  of  depth  300m 
with  c  =  1500m/s,  we  place  two  50Hz  sources  in  phase  200m  apart  both  at  z  —  25m.  The  point 


3DPE  (solid  curves)  and  with  the  2DPE  (dashed  curves)  for  (a)  6  =  0°,  (b)  8  =  20°,  (c)  6  =  40°, 


(d)  6  =  60°,  (e)  8  =  80°,  (Figure  2.  cont’d  on  next  page) 


Loss  (dB  re  lm)  Loss  (dB  re  1m) 
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Thela  =  180  (leg  The  la  =  1G0  deg 


Thela  =  140  deg 


Range  (km) 


Thela  =  hdO  deg 


Thela  =  100  deg 


Figure  2  cont’d:  Transmission  loss  at  z  =  30m  in  an  ocean  with  a  corrugated  bottom  generated 
with  the  3DPE  (solid  curves)  and  with  the  2DPE  (dashed  curves)  for  (f)  0  =  100°,  (g)  0  =  120°, 
(h)  9  =  140°,  (i)  9  =  160°,  (j)  9  =  180°. 
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midway  between  the  sources  is  r  =  0.  In  the  sediment,  c  =  1600m/s,  p  =  1.5g/cm,  and  /3  =  0.5. 
The  homogeneous  half-space  field  [12] 

Ph  =  -7-  exp(iM-)  -  -T-  exp  (ikad+)  (2.10) 

«+ 

=  r2  +  (z±z0)2  (2.11) 

is  used  as  an  initial  condition  at  r  =  150m,  and  the  3DPE  is  applied  to  march  the  field  out  to  r 
=  300m.  The  modal  radiation  patterns  [13]  201og]0  |<  P,\j> j  >|  and  20logi0  j<  Ph,il>i  >|  appear 
in  Figure  1,  where  V'i  is  the  first  waveguide  normal  mode  and  <,>  is  the  inner  product  associated 
with  the  depth  operator  of  Eq.  (2.1).  The  half-space  radiation  pattern  evolves  substantially  over 
150m  <  r  <  300m,  and  the  radiation  patterns  are  in  excellent  agreement  at  r  =  300m. 

To  demonstrate  that  azimuthal  diffraction  can  be  important,  we  consider  a  25Hz  source  at  z  = 
25m  in  an  ocean  of  depth 


d  = 


3 


(2.12) 


which  depends  only  on  the  Cartesian  coordinate  x  =  rcos0.  The  maximum  slope  of  the  ocean 
bottom  is  about  3°.  In  the  water  column,  we  take  c  —  1500m/s.  In  the  sediment,  c  =  1700m/s,  p 
=  1.5g/cm3,  and  fi  =  0.5.  The  domain  is  truncated  at  z  =  600m,  and  the  attenuation  is  increased 
linearly  to  20  in  the  lower  100m  of  the  domain  to  prevent  reflections.  The  grid  spacings  are  Ar 
=  5m,  Az  =  lm,  and  A 9  =  0.25°.  Transmission  loss  computed  with  the  2PDE  and  the  3DPE 
appears  in  Figure  2.  We  observe  that  the  solutions  are  in  good  agreement  for  9  near  0°  and  180°. 
Near  90°,  however,  there  is  a  large  difference  between  the  solutions  because  energy  gets  trapped  in 
the  deep  parts  of  the  corrugated  bottom  and  channeled  in  the  y-direction.  Similar  effects  should 
be  important  for  many  realistic  shallow-water  problems. 

The  three-dimensional  behavior  we  have  observed  can  also  be  illustrated  with  ray  tracing.  The 
grazing  angle  <j>  of  a  ray  is  defined  to  be  the  angle  the  ray  makes  with  the  ocean  surface.  Rays  are 
traced  from  the  origin  for  <j>  =  2n°  for  1  <  n  <  15.  Rays  that  encounter  the  ocean  bottom  reflect 
(terminate)  if  incident  within  (beyond)  the  critical  angle.  Projections  of  the  ray  paths  onto  the 
ocean  surface  appear  in  Figure  3.  Ray  channeling  is  significant  only  for  60°  <  9  <  120°,  which  is 
consistent  with  the  3DPE  results.  For  the  case  of  9  =  80°,  the  rays  in  Figure  3  are  trapped  in  the 
channel  for  5  <  n  <  13. 


3.  Rough  surface  modeling 


In  the  farfield,  we  consider  the  plane  wave 


Pi  =  exp  [ik(rcos<j>  -  zsin<£)] 


(3.1) 


incident  upon  a  rough  ocean  surface.  As  in  Ref.  7,  we  assume  that  the  scattered  field  is  given  bv 
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P,  =  -R(<t>)  exp  [ifc(r  cos<£  +  zsin  <£)] ,  (3.2) 

where  R  is  the  reflection  coefficient.  Since  a.  rough  ocean  surface  results  in  variations  in  path  lengths 
and  thus  phase  distortions,  we  allow  R  to  be  complex.  Since  a  rough  surface  behaves  like  a  smooth 
reflector  at  small  grazing  angles,  [14]  fi(0)  =  1.  Since  |  < j>  |<<  1  for  rays  in  the  farfield,  we  may 
assume 


R  ~  1  +  +  —  (3.3) 

We  use  Eqs.  (3.1),  (3.2),  and  (3.3)  to  determine  (3j  such  that  the  following  expression  for  the  total 
field  P  =  +  P,  is  correct  to  the  appropriate  order  in  <f>  at  z  =  0: 


A  P  +  h  +  A  I?  +  Tr  (A  P  +  A  +  A  H)  =0.  (3.4) 


To  avoid  a  degenerate  solution,  it  is  necessary  to  assume  that  some  of  the  coefficients  in  Eq.  (3.4 ) 
vanish.  For  the  case  0j  =  0  for  j  >  2,  we  obtain 


2kP  +  iai-~ -  =0.  (3.5) 

OZ 

This  approach  can  also  be  applied  for  the  cases  /?,  =  =  /J6  =  0  (wide  angle)  and  /J,  =  (i2  =  0 

(very  wide  angle)  for  better  accuracy.  To  model  rough  surfaces  with  the  3DPE,  we  replace  the 
pressure  release  boundary  condition  P  =  0  that  is  usually  used  in  PE  models  with  Eq.  (3.4). 
This  approach  is  easily  implemented  numerically  into  Eq.  (2.7).  A  nonphysical  depth  grid  point  is 
introduced  at  z  =  —A z  and  Eq.  (3.4)  is  assigned  to  this  point.  Equation  (3.4)  is  discretized  using 
centered  differences  in  both  r  and  z,  and  Eq.  (2.7)  is  assigned  to  each  of  the  physical  depth  grid 
points  z  >  0  and  discretized  with  Galerkin’s  method  and  Crank  Nicolson  integration  as  described 
in  Ref.  10. 

The  plane  wave  reflection  coefficient  can  be  approximately  introduced  into  the  normal  mode 
solution  [15] 


p  =  5Z  aj^i(2)exP  (ikjr) 
j 


(3.6) 
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liaiigo  (kin) 

Figure  4:  Transmission  loss  at  z  =  25m  in  an  ocean  with  a  rough  surface.  The  PE  solution  (solid 
curve)  versus  the  normal  mode  solution  (dashed  curve).  The  broken  curve  is  the  smooth  surface 
PE  solution. 

as  a  perturbation,  where  aj  are  constants,  if>j  are  modes,  and  kj  are  complex  eigenvalues.  We  assume 
that  the  leading-order  effect  of  this  perturbation  is  to  increase  loss.  For  simplicity,  c  is  taken  to 
be  constant  in  the  water  column.  The  number  of  surface  bounces  that  a  mode  propagating  at  the 
angle  <j>j  experiences  over  the  range  r  in  an  ocean  of  depth  d  is  rtan<£y/2d.  Thus  a  leading-order 
rough  surface  normal  mode  solution  is 


P  =  Y,  a^(z)R(4>j)  rlan^/2dexp  (ikjr).  (3.7) 

3 

Equation  (3.7)  is  similar  to  the  rough  surface  normal  mode  model  of  Ref.  8,  which  was  used  to 
model  both  rough  ocean  surfaces  and  rough  ocean  bottoms.  Thus  the  rough  surface  3DPE  should 
be  valid  even  for  modeling  rough  ocean  bottoms.  Since  the  boundary  condition  is  applied  at  the 
ocean  surface,  this  is  not  obvious  from  the  derivation. 

To  test  the  rough  surface  PE,  we  consider  a  stratified  ocean  of  depth  200m  in  which  c  =  1500m/s. 
We  take  at  =  -1  and  place  a  50Hz  source  at  z  =  25m.  In  the  ocean  bottom,  c  =  1600m/s,  p  = 
1.5g/cm3,  and  0  =  0.5.  Transmission  loss  computed  with  Eq.  (3.7)  and  with  the  rough  and  smooth 
surface  2DPE  models  appears  in  Figure  4.  The  rough  surface  solution  exhibits  more  loss  than  the 
smooth  surface  solution.  Furthermore,  the  rough  surface  2DPE  and  normal  mode  solutions  are 
in  good  agreement.  The  agreement  is  not  perfect,  however,  because  the  normal  mode  model  was 
derived  from  the  PE  model  by  approximation.  Since  Oj  was  taken  to  be  real,  the  minima  of  the 
transmission  loss  curves  are  coincident  for  the  smooth  surface  and  rough  surface  results. 

To  illustrate  the  rough  surface  PE  model  in  range-dependent  domains,  we  consider  the  cylindri- 
cally  symmetric  rough  bottom  problem  of  Ref.  6.  The  ocean  bottom  is  smooth  for  r  <  5km.  In  the 
rough  bottom  region  5km  <  r  <  rjif,d(r)  is  a  square  wave  with  d  =  90m  at  the  50m  wide  minima 
and  d  =  100m  at  the  50m  wide  maxima.  The  cases  r\ i  =  10km  and  r\j  =  o o  are  considered.  A 
25Hz  source  is  placed  at  z  =  18m.  In  the  bottom,  c  =  1704.5m/s,  p  =  2.5g/cm3,  and  0  =  0.5.  It 
is  clear  from  Figures  4  and  5  of  Ref.  6  that  the  rough  bottom  alters  both  phases  and  amplitude  as 
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Figure  5:  Transmission  loss  at  z  =  50m  in  an  ocean  with  a  rough  surface  for  (a)  r  >  5km  and  (b) 
for  5km  <  r  <  10km.  The  solid  curve  is  the  rough  surface  PE  solution.  The  dashed  curve  is  the 
smooth  surface  PE  solution. 

the  minima  of  the  transmission  loss  curves  for  the  smooth  and  rough  cases  do  not  coincide.  Thus 
R  is  complex  for  this  problem.  Transmission  loss  for  the  smooth  surface  PE  solution  and  the  rough 
surface  PE  solution  for  at  =  -|exp(2xt/9)  appear  in  Figure  5.  The  discontinuities  in  O]  were 
handled  with  linear  variations  over  short  ranges.  The  rough  surface  PE  solutions  are  similar  to  the 
normal  mode  solutions  in  Ref.  6. 


4.  Efficient  PE  algorithms 


Finite  difference  solutions  of  a2DPE  involve  a  system  of  linear  equations  involving  a  tridiagonal 
matrix.  It  is  natural  to  solve  this  system  using  Gaussian  elimination  from  the  ocean  surface  down 
to  the  bottom  of  the  grid.  Since  the  system  is  solved  repeatedly  as  the  solution  is  marched  in  range, 
a  great  deal  of  efficiency  can  be  gained  by  decomposing  the  matrix  into  upper  and  lower  triangular 
matrices.  An  efficient  2DPE  code  for  range-independent  problems  would  have  FORTRAN  loops  in 
the  tridiagonal  solver  subroutine  similar  to: 


DO  1  I  =  2,  N 

U(I)  =  A(I)*U(I)  +  B(I)*U(1-1) 

1  CONTINUE 

DO  2  I  =  N-l,  1,  -1 

U(I)  =  U(I)  +  C(I)*U(I+1) 

2  CONTINUE 
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The  first  loop  corresponds  to  downward  elimination.  The  second  loop  corresponds  to  back  substi¬ 
tution. 

The  diagonals  A,  B,  and  C  are  determined  from  the  triangular  decomposition  and  stored. 
Since  computers  perform  multiplication  significantly  faster  than  division,  it  is  better  to  store  these 
constants  as  factors  rather  than  divisors  as  we  illustrate  with  a  benchmark  problem.  [16,17]  We 
consider  a  range-dependent  problem  for  which  the  ocean  depth  decreases  linearly  from  200m  at  r 
=  0  to  zero  at  r  =  4km.  A  25IIz  source  is  placed  at  z  =  100m  in  an  ocean  in  which  c  =  1500m/s. 
In  the  bottom,  c  =  1700m/s,  p  =  1.5g/cm3,  and  /3  =0.5.  We  truncate  the  domain  at  z  =  2km 
and  take  A r  =  5m  and  Az  =  0.5m.  These  parameters  were  used  in  Ref.  16  to  solve  this  problem 
using  the  PE  code  IFD,  [18,19]  which  performs  division  in  the  loops.  The  run  time  reported  in  Ref. 
16  was  7min  on  an  FPS-164  computer.  Using  Gaussian  elimination  with  the  PE  code  FEPE,  [20] 
which  does  not  have  division  in  the  solver  loops,  the  run  time  is  2min  on  an  FPS-164. 

Gaussian  elimination  is  not  the  most  efficient  algorithm  for  problems  involving  range-dependent 
ocean  depth.  As  depth  varies,  it  is  necessary  to  modify  the  entries  of  the  matrix  rows  corresponding 
to  the  vicinity  of  the  ocean  bottom  interface.  Thus  it  is  necessary  to  repeat  downward  elimination 
below  the  interface.  A  more  efficient  approach  is  to  eliminate  the  entries  below  the  main  diagonal 
from  the  ocean  surface  down  to  the  ocean  bottom  and  to  eliminate  the  entries  above  the  main 
diagonal  from  the  bottom  of  the  domain  up  to  the  ocean  bottom.  With  this  approach,  it  is 
necessary  to  repeat  the  elimination  process  only  for  the  rows  near  the  ocean  bottom  as  depth  varies, 
and  run  time  is  essentially  independent  of  bathymetry  variations.  Furthermore,  this  algorithm  is 
vectorizable  and  may  improve  run  times  by  up  to  a  factor  of  two  on  a  vector  machine.  Using  this 
approach  in  FEPE,  the  run  time  for  the  benchmark  problem  is  lmin  on  an  FPS-164.  For  rough 
surfaces  in  which  the  require  frequent  updates,  it  may  be  useful  to  perform  Gaussian  elimination 
from  the  bottom  of  the  grid  up  to  the  ocean  surface. 

5.  Conclusion 

A  3DPE  that  handles  wide  angles  in  depth  and  narrow  angels  in  azimuth  and  rough  surfaces  has 
been  derived  and  solved  numerically  using  alternating  directions.  Surface  roughness  was  handled 
by  assuming  a  reflection  coefficient  that  depends  on  grazing  angle  and  implemented  into  the  3DPE 
as  a  homogeneous  boundary  condition.  The  accuracy  of  the  3DPE  was  tested  with  a  benchmark 
calculation.  Calculations  were  presented  to  illustrate  the  rough  surface  model  and  to  demonstrate 
that  azimuthal  diffraction  can  be  important  in  shallow  water.  Efficient  PE  algorithms  have  been 
considered. 
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